NewGridIntegerFromNetCDF Subroutine

private subroutine NewGridIntegerFromNetCDF(layer, fileName, variable, stdName, time)

Uses

create a new grid_integer reading data from NetCDF file The variable to read can be defined by its current name or the standard_name. The dimensions x and j of the variable is calculated searching from the dimensions of the couple of variables with 'standard_name' equal to 'projection_x_coordinate' and 'projection_y_coordinate' for projected reference systems or 'longitude' and 'latitude' for geographic reference systems or 'grid_longitude' and 'grid_latitude' for rotated pole systems If a comprehensible reference systems is not found a geodetic reference system is supposed. Once the variable is retrieved, offset and scale factor are applied and a check on minimum and maximum valid value is performed.

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(out) :: layer

gridreal to return
character(len=*), intent(in) :: fileName

NetCDF file to read
character(len=*), intent(in), optional :: variable

variable to read
character(len=*), intent(in), optional :: stdName

standard name of the variable to read
type(DateTime), intent(in), optional :: time

time of the grid to read

Variables

Type Visibility Attributes Name Initial
character(len=80), public :: attribute
integer(kind=short), public :: i

loop index
integer(kind=short), public :: ios

error return code
integer(kind=short), public :: j

loop index
integer(kind=short), public :: nDimsVar

number of dimensions of a variable
integer(kind=short), public :: nVars

number of variables
integer(kind=short), public :: ncId

NetCdf Id for the file
integer(kind=short), public :: ncStatus

error code return by NetCDF routines
integer(kind=long), public :: offset
integer(kind=long), public :: scale_factor
character(len=1), public :: shp
integer, public, ALLOCATABLE :: slice(:)
integer(kind=long), public, ALLOCATABLE :: tempGrid(:,:)
integer(kind=short), public :: varId

variable Id
character(len=100), public :: variableName

Source Code

SUBROUTINE NewGridIntegerFromNetCDF &
!
(layer, fileName, variable, stdName, time)

USE StringManipulation, ONLY: &
!Imported routines:
StringCompact


IMPLICIT NONE

! Scalar arguments with intent(in):
CHARACTER (LEN = *), INTENT(in) :: fileName  !!NetCDF file to read
CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: variable  !!variable  to read
CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: stdName  !!standard name of 
                                                      !!the variable  to read
TYPE (DateTime), OPTIONAL, INTENT(in) :: time  !!time of the grid to read

! Arguments with intent(out):
TYPE (grid_integer), INTENT (out)          :: layer  !!gridreal to return

! Local scalars:
INTEGER (KIND = short)          :: ios !!error return code
INTEGER (KIND = short)          :: ncStatus !!error code return by NetCDF routines
INTEGER (KIND = short)          :: ncId  !!NetCdf Id for the file
INTEGER (KIND = short)          :: varId !!variable Id
INTEGER (KIND = short)          :: nVars !!number of variables
CHARACTER (LEN = 80)            :: attribute
INTEGER (KIND = short)          :: nDimsVar !!number of dimensions of a variable
INTEGER (KIND = short)          :: i, j   !!loop index
CHARACTER (LEN = 100)           :: variableName
CHARACTER (LEN = 1)             :: shp
INTEGER (KIND = long)           :: scale_factor
INTEGER (KIND = long)           :: offset


! Local arrays:
INTEGER, ALLOCATABLE            :: slice (:)
INTEGER (KIND = long), ALLOCATABLE  :: tempGrid (:,:) 

!------------end of declaration------------------------------------------------

!------------------------------------------------------------------------------
![1.0] open NetCDF dataset with read-only access
!------------------------------------------------------------------------------
ncStatus = nf90_open (fileName, NF90_NOWRITE, ncId)
IF (ncStatus /= nf90_noerr) THEN
  CALL Catch ('error', 'GridLib',        &
  TRIM (nf90_strerror (ncStatus) )//': ',  &
  code = ncIOError, argument = fileName )
ENDIF

!------------------------------------------------------------------------------
![2.0] get x and y size and allocate array
!------------------------------------------------------------------------------

CALL GetXYSizesFromFile (ncId, layer % jdim, layer % idim, shape = shp)

!allocate grid
ALLOCATE ( layer % mat (layer % idim, layer % jdim), STAT = ios )
IF (ios /= 0) THEN
  CALL Catch ('error', 'GridLib',  &
  'memory allocation ',  &
  code = memAllocError,argument = variable )
ENDIF

!allocate temporary grid
ALLOCATE ( tempGrid (layer % jdim, layer % idim), STAT = ios )
IF (ios /= 0) THEN
  CALL Catch ('error', 'GridLib',  &
  'memory allocation ',  &
  code = memAllocError,argument = variable )
ENDIF

!------------------------------------------------------------------------------
![3.0] get values
!------------------------------------------------------------------------------
IF ( PRESENT (variable) ) THEN
  ncStatus = nf90_inq_varid (ncId, variable, varId)
  CALL ncErrorHandler (ncStatus)
ELSE IF (PRESENT (stdName) ) THEN !search variable corresponding to standard name

  !inquire dataset to retrieve number of dimensions, variables 
  !and global attributes
  ncStatus = nf90_inquire(ncId,  nVariables = nVars )
                          
  CALL ncErrorHandler (ncStatus)
  
  DO i = 1, nVars
    attribute = ''
    ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', &
                             values = attribute)
    IF (ncStatus == nf90_noerr) THEN !error is raised if 'standard_name' is not found
      IF ( StringCompact(attribute(1:LEN_TRIM(attribute))) == stdName(1:LEN_TRIM(stdName)) )THEN
         varId = i
      END IF
    END IF
  END DO
END IF

!verify number of dimensions of variable to read
ncstatus = nf90_inquire_variable(ncId, varId, ndims = nDimsVar) 
CALL ncErrorHandler (ncStatus)

!If number of dimensions = 3 read info about time
IF (nDimsVar == 3) THEN
  ALLOCATE ( slice(3) )
  !Read time informations
  CALL ParseTime (ncId, layer % reference_time, layer % time_unit)

  !Define current time slice (if present)
  IF (PRESENT (time) ) THEN
    slice = (/ 1, 1, 1 /)
    slice(3) = TimeIndex (ncId, layer % reference_time, &
                          layer % time_unit, time)
    layer % current_time = time
    layer % time_index = slice(3)
  ELSE
    slice = (/ 1, 1, 1 /)
  ENDIF
ELSE
   ALLOCATE ( slice(2) )
   slice = (/ 1, 1/)
END IF 
ncStatus = nf90_get_var (ncId, varId, tempGrid , start = slice)
CALL ncErrorHandler (ncStatus)

!transpose temporary matrix to grid specification
CALL SwapGridIntegerForward (tempGrid, layer % mat)

!deallocate temporary grid
DEALLOCATE (tempGrid)

!------------------------------------------------------------------------------
![4.0] get attributes
!------------------------------------------------------------------------------

ncStatus = nf90_get_att (ncId, varId, name = 'standard_name', &
                             values = layer % standard_name)
                            
ncStatus = nf90_get_att (ncId, varId, name = 'long_name', &
                             values = layer % long_name)                              
                                                       
ncStatus = nf90_get_att (ncId, varId, name = 'units', &
                             values = layer % units) 
                            
ncStatus = nf90_get_att (ncId, varId, name = 'varying_mode', &
                             values = layer % varying_mode)
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) layer % varying_mode = 'sequence'                            
                               
ncStatus = nf90_get_att (ncId, varId, name = '_FillValue', &
                             values = layer % nodata)

!if _FillValue is not defined search for missing_value   
IF (ncStatus /= nf90_noerr)  THEN                         
  ncStatus = nf90_get_att (ncId, varId, name = 'missing_value', &
                             values = layer % nodata)
END IF   

 !if missing_value is not defined set to default
IF (ncStatus /= nf90_noerr)  THEN                         
  layer % nodata = MISSING_DEF_REAL
END IF   
                             
ncStatus = nf90_get_att (ncId, varId, name = 'valid_min', &
                             values = layer % valid_min) 
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) layer % valid_min = layer % nodata 
                                                         
ncStatus = nf90_get_att (ncId, varId, name = 'valid_max', &
                             values = layer % valid_max) 
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) layer % valid_max = layer % nodata
                                                        
ncStatus = nf90_get_att (ncId, varId, name = 'scale_factor', &
                             values = scale_factor)
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) scale_factor = 1
                             
ncStatus = nf90_get_att (ncId, varId, name = 'add_offset', &
                             values = offset)
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) offset = 0  

ncStatus = nf90_get_att (ncId, varId, name = 'esri_pe_string', &
                             values = layer % esri_pe_string)
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) layer % esri_pe_string = ''                         

!------------------------------------------------------------------------------
![5.0] check values and apply scale factor and offset
!------------------------------------------------------------------------------
!retrieve name of variable
ncstatus = nf90_inquire_variable(ncId, varId, name = variableName) 
layer % var_name = variableName
                              
DO i = 1, layer % jdim
  DO j = 1, layer % idim 
    IF ( layer % mat (i,j) /= layer % nodata ) THEN
    
      !apply scale factor
      layer % mat (i,j) = layer % mat (i,j) * scale_factor
      
      !add offset
      layer % mat (i,j) = layer % mat (i,j) + offset
      
      !check lower bound
      IF ( layer % valid_min /= layer % nodata .AND. &
           layer % mat (i,j) < layer % valid_min ) THEN
        layer % mat (i,j) = layer % valid_min
        CALL Catch ('info', 'GridLib', 'corrected value exceeding &
         lower bound in variable: ', argument = variableName )
      END IF  
    
      !check upper bound
      IF ( layer % valid_max /= layer % nodata .AND. &
           layer % mat (i,j) > layer % valid_max ) THEN
        layer % mat (i,j) = layer % valid_max
        CALL Catch ('info', 'GridLib', 'corrected value exceeding &
         upper bound in variable: ', argument = variableName )
      END IF  
    END IF
  END DO
END DO   

!------------------------------------------------------------------------------
![6.0] set file name
!------------------------------------------------------------------------------
layer % file_name = fileName

!------------------------------------------------------------------------------
![7.0] read georeferencing informations from netCDF file
!------------------------------------------------------------------------------

IF (shp == 'm') THEN !2 dimensional (matrix) coordinate for not regular grid
  IF (.NOT. ALLOCATED (layer % Iraster)) THEN
    ALLOCATE (layer % Iraster(layer%idim,layer%jdim))
    ALLOCATE (layer % Jraster(layer%idim,layer%jdim))
    CALL GetGeoreferenceFromNCdataSet (ncId, varId, layer % cellsize, &
                 layer % xllcorner, layer % yllcorner, &
                 layer % grid_mapping, layer % Iraster, layer % Jraster, .TRUE.)
  END IF
  
  
  !create a temporary copy of layer % mat
  IF (ALLOCATED (tempGrid) ) THEN
    DEALLOCATE (tempGrid)
    ALLOCATE (tempGrid(layer%idim,layer%jdim))
  ELSE
    ALLOCATE (tempGrid(layer%idim,layer%jdim))
  END IF
  tempGrid = layer % mat
  layer % mat = layer % nodata
  !fill in the regular grid
  DO i = 1, layer%idim
    DO j = 1, layer%jdim
      IF (layer % Iraster (i,j) /= -9999 .AND.  layer % Jraster (i,j) /= -9999) THEN
        layer % mat (layer%Iraster(i,j),layer%Jraster(i,j)) = tempGrid (i,j)
      ELSE
        layer % mat(i,j) = layer % nodata
      END IF
    END DO
  END DO
  DEALLOCATE (tempGrid)
ELSE !regular grid stores coordinate in 1dimensional array (vector)
  CALL GetGeoreferenceFromNCdataSet (ncId, varId, layer % cellsize, &
                 layer % xllcorner, layer % yllcorner, &
                 layer % grid_mapping, layer % Iraster, layer % Iraster, .FALSE.)                                
END IF
                                                                                                                             
!------------------------------------------------------------------------------
![8.0] close NetCDF dataset
!------------------------------------------------------------------------------
ncStatus = nf90_close (ncid)
CALL ncErrorHandler (ncStatus)

END SUBROUTINE NewGridIntegerFromNetCDF